program gen_mesh

    use mesh
    use mesh_file
    
    implicit none
    
    integer i, j
    character(100) file_name
    
    is_domain_regular = .true.
    is_mesh_uniform = .true.
    
    write(*, *) "************ GMIT Regular Mesh Generator ************"
    write(*, *) "Available mesh class:"
    write(*, *) "1. Uniform line mesh (1D)"
    write(*, *) "2. Uniform quadrilateral mesh (2D)"
    write(*, *)
    write(*, *) "Available boundary class:"
    write(*, *) "1. Periodic boundary"
    write(*, *)
    write(*, *) "Input mesh class:"
    write(*, "('> ')", advance="no")
    read(*, *) mesh_class
    
    select case(mesh_class)
    case(1)
        num_dim = 1
    case(2)
        num_dim = 2
    end select
    allocate(num_grid_(num_dim))
    allocate(num_cell_(num_dim))
    allocate(boundary_x(2,num_dim))
    allocate(span(num_dim))
    allocate(boundary_class(2,num_dim))
    allocate(dx(num_dim))
    
    write(*, *) "Input number of grids:"
    write(*, "('> ')", advance="no")
    read(*, *) num_grid_
    
    write(*, *) "Input boundary coordinates:"
    write(*, "('> ')", advance="no")
    read(*, *) boundary_x
    
    write(*, *) "Input boundary classes:"
    write(*, "('> ')", advance="no")
    read(*, *) boundary_class
    
    write(*, *) "Input output mesh file name:"
    write(*, "('> ')", advance="no")
    read(*, *) file_name
    
    ! preprocess input
    do i = 1, num_dim
        span(i) = boundary_x(2,i)-boundary_x(1,i)
        if(boundary_class(1,i) == 1 .and. boundary_class(2,i) == 1) then
            dx(i) = span(i)/num_grid_(i)
        end if
    end do
    
    select case(mesh_class)
    case(1) ! line
        ! Grid
        ! memory allocation
        num_grid = num_grid_(1)
        allocate(grids(num_grid))
        do i = 1, num_grid
            allocate(grids(i)%x(num_dim))
        end do
        ! setting of grid coordinate
        if(boundary_class(1,1) == 1 .and. boundary_class(2,1) == 1) then
            do i = 1, num_grid
                grids(i)%x(1) = i*dx(1)
            end do
        end if
        ! Internal Cell
        if(boundary_class(1,1) == 1 .and. boundary_class(2,1) == 1) then
            ! memory allocation
            num_cell_ = num_grid_
            num_cell = num_grid
            allocate(cells(num_cell))
            do i = 1, num_cell
                cells(i)%num_vertex = 2
                allocate(cells(i)%vertex_index(cells(i)%num_vertex))
                allocate(cells(i)%vertex_x(num_dim,cells(i)%num_vertex))
                cells(i)%num_neighb_cell = 2
                allocate(cells(i)%neighb_cell_index(cells(i)%num_neighb_cell))
                allocate(cells(i)%is_neighb_cell_far( &
                    num_dim,cells(i)%num_neighb_cell))
                cells(i)%num_edge = 2
                allocate(cells(i)%edges(cells(i)%num_edge))
                do j = 1, cells(i)%num_edge
                    cells(i)%edges(j)%num_vertex = 1
                    allocate(cells(i)%edges(j)%local_vertex_index( &
                        cells(i)%edges(j)%num_vertex))
                    allocate(cells(i)%edges(j)%norm(num_dim))
                end do
            end do
            ! set grid connectivity
            do i = 1, num_cell
                cells(i)%vertex_index(1) = i-1
                cells(i)%vertex_index(2) = i
            end do
            cells(1)%vertex_index(1) = num_cell
            ! solidify grid coordinates avoiding boundary effect
            do i = 1, num_cell
                cells(i)%vertex_x(:,1) = grids(cells(i)%vertex_index(1))%x
                cells(i)%vertex_x(:,2) = grids(cells(i)%vertex_index(2))%x
            end do
            cells(1)%vertex_x(1,1)= span(1)-grids(cells(1)%vertex_index(1))%x(1)
            ! set cell connectivity
            do i = 1, num_cell
                cells(i)%neighb_cell_index(1) = i-1
                cells(i)%is_neighb_cell_far(1,1) = .false.
                cells(i)%neighb_cell_index(2) = i+1
                cells(i)%is_neighb_cell_far(1,2) = .false.
            end do
            cells(1)%neighb_cell_index(1) = num_cell
            cells(1)%is_neighb_cell_far(1,1) = .true.
            cells(num_cell)%neighb_cell_index(2) = 1
            cells(num_cell)%is_neighb_cell_far(1,2) = .true.
            ! set edge
            do i = 1, num_cell
                cells(i)%edges(1)%local_vertex_index(1) = 1
                cells(i)%edges(1)%norm = [-1.0d0]
                cells(i)%edges(2)%local_vertex_index(1) = 2
                cells(i)%edges(2)%norm = [1.0d0]
            end do
        end if
        ! Boundary Cell (no need for regular domain)
    case(2) ! quadrilateral
        ! Grid
        
        ! Cell
    end select

    call mesh_file_writer(file_name)

    stop

end program gen_mesh